Development of a Monte Carlo-wave model to simulate time domain diffuse correlation spectroscopy measurements from first principles

Abstract. Significance Diffuse correlation spectroscopy (DCS) is an optical technique that measures blood flow non-invasively and continuously. The time-domain (TD) variant of DCS, namely, TD-DCS has demonstrated a potential to improve brain depth sensitivity and to distinguish superficial from deeper blood flow by utilizing pulsed laser sources and a gating strategy to select photons with different pathlengths within the scattering tissue using a single source–detector separation. A quantitative tool to predict the performance of TD-DCS that can be compared with traditional continuous wave DCS (CW-DCS) currently does not exist but is crucial to provide guidance for the continued development and application of these DCS systems. Aims We aim to establish a model to simulate TD-DCS measurements from first principles, which enables analysis of the impact of measurement noise that can be utilized to quantify the performance for any particular TD-DCS system and measurement geometry. Approach We have integrated the Monte Carlo simulation describing photon scattering in biological tissue with the wave model that calculates the speckle intensity fluctuations due to tissue dynamics to simulate TD-DCS measurements from first principles. Results Our model is capable of simulating photon counts received at the detector as a function of time for both CW-DCS and TD-DCS measurements. The effects of the laser coherence, instrument response function, detector gate delay, gate width, intrinsic noise arising from speckle statistics, and shot noise are incorporated in the model. We have demonstrated the ability of our model to simulate TD-DCS measurements under different conditions, and the use of our model to compare the performance of TD-DCS and CW-DCS under a few typical measurement conditions. Conclusion We have established a Monte Carlo-Wave model that is capable of simulating CW-DCS and TD-DCS measurements from first principles. In our exploration of the parameter space, we could not find realistic measurement conditions under which TD-DCS outperformed CW-DCS. However, the parameter space for the optimization of the contrast to noise ratio of TD-DCS is large and complex, so our results do not imply that TD-DCS cannot indeed outperform CW-DCS under different conditions. We made our code available publicly for others in the field to find use cases favorable to TD-DCS. TD-DCS also provides a promising way to measure deep brain tissue dynamics using a short source–detector separation, which will benefit the development of technologies including high density DCS systems and image reconstruction using a limited number of source–detector pairs.


Introduction
Cerebral blood flow (CBF) is an important indicator of brain function and health. [1][2][3][4][5][6] Diffuse correlation spectroscopy (DCS) serves as a major optical technique that is capable of measuring blood flow non-invasively and continuously. [7][8][9][10][11] In DCS, coherent light is incident on the surface of the scattering medium, and the re-emitted scattered light is collected by a detector at a certain distance away from the source position, typically in the range of 10 to 30 mm. Dynamics in the tissue, mainly arising from blood flow induces phase variations of the scattered light that alters the interference pattern of the partial waves of the re-emitted light, namely, the speckle pattern, thus causing light intensity fluctuations with time at the detector. It is common to quantify these temporal fluctuations using a temporal intensity autocorrelation function. The decay rate of this autocorrelation function provides a measure of the blood flow index. Numerical and theoretical tools have been developed to guide interpretations of experimental measurements. For traditional continuous-wave DCS (CW-DCS), the analytical expressions of the autocorrelation functions have been well established in diffusion theory, 7,12,13 and Monte Carlo simulations have been developed to numerically obtain the field autocorrelation function for inhomogeneous media. 8,14 Together with the noise model obtained analytically, 15 the performance of a CW-DCS system for a particular measurement geometry can be theoretically predicted.
Recently, a time-domain variant of DCS (TD-DCS) has been developed, which utilizes a pulsed laser and a gating strategy to select photons detected at different arrival times. [16][17][18][19][20] This technology has the potential to differentiate blood flow information from different layers, such as the skull and the brain tissue using a single source-detector separation. The theoretical basis for TD-DCS has been established taking into account the effects of the instrument response function (IRF), which includes the incident pulse shape and the broadening of the pulse arising from the instrument, and the coherence length of the light source on the shape of the intensity autocorrelation function. 18 Ideally, a narrower width of the IRF will result in better specificity to different photon pathlengths. However, the width of IRF is fundamentally limited by the coherence length, and a narrow pulse in time will result in a short coherence length, which will degrade signal to noise ratio (SNR) in TD-DCS measurements. 18 A pioneering study on the performance of TD-DCS has been conducted assuming the same noise model 15 that is developed for CW-DCS. 21 A proper noise model that quantifies the performance of TD-DCS measurements in terms of SNR and contrast to noise ratio (CNR) for any IRF, coherence length, measurement geometry, and brain activation pattern is still lacking.
We have constructed a numerical model which is capable of simulating both CW-and TD-DCS measurements from first principles. We have integrated Monte Carlo simulations with wave calculations to obtain the intensity fluctuations induced by tissue dynamics. The speckle intensity is converted to photon counts within each bin time mimicking the experimental measurements. The numerical results are compared with existing theoretical expressions of the average intensity autocorrelation function for CW-DCS 14 and TD-DCS. 18 The simulated noise is also validated with the existing noise model for CW-DCS. 15 We have demonstrated the value of using our model to predict the performance for any TD-DCS systems and measurement geometry in terms of CNR due to brain activation. We have found that for a given set of parameters, brain geometry, and the IRF profile utilized in this paper, CW-DCS at a large source-detector separation (30 mm) out-performs TD-DCS. We can increase the detected photon flux per bin time of TD-DCS at a small source-detector separation (10 mm) by increasing the gate width at the cost of reduced specifity to a particular photon pathlength or depth to achieve a performance comparable to CW-DCS at a medium source-detector separation (20 mm). We could not find realistic measurement conditions under which TD-DCS outperformed CW-DCS. However, since the performance of TD-DCS is highly dependent on the IRF profile 21 and the other measurement parameters, the quantitative results predicted in this paper can be different for researchers using different input parameters in our model. The performance for a given TD-DCS system can be easily calculated by changing the parameters of the numerical model, which is made publicly available. 22 TD-DCS is also capable of differentiating blood flow at different depths using a single source-detector separation. This model provides guidance for the experimental development of novel TD-DCS systems.

Theory
DCS quantifies the dynamics of the sample by calculating the autocorrelation function of the speckle intensity fluctuations of light reemitted from tissue. For a single scattering event shown in Fig. 1(a), the electric field measured at the detector is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 2 6 4Ẽ Here,Ẽ is the electric field reaching the detector located atR d ;r n is the position of the n'th scatterer; N is the total number of the scatterers;Ẽ 0 is the incident field; andq n ¼k out −k in is the momentum transfer, wherek out andk in are the output and incident wave vectors, respectively. We only consider elastic scattering such that jk out j ¼ jk in j ¼ k.
For the case of multiple light scattering, Eq. (1) is revised as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 1 1 2Ẽ expð−iq ns ·r ns Þ: The geometry of the Monte Carlo simulation used in this manuscript and the illustration of the gating strategy for TD-DCS. The size of the detector is 2 mm, the thickness of the first layer or first tissue type is 15 mm, the reduced scattering coefficient is μ s 0 ¼ 1 mm −1 and absorption coefficient μ a ¼ 0.01 mm −1 . The optical properties μ s 0 and μ a are set to be the same for the two layers. The dynamics in the second layer is increased when brain activation is introduced. The source detector separation is set to be ρ ¼ 20 mm unless otherwise stated. Measurements at a later gate will be able to select photons with longer pathlengths, thus it is more sensitive to dynamics at deeper depths. Here, s denotes the s'th scattering event, and N sn is the total number of scattering events for the n'th component of the wave, or the n'th photon in our Monte Carlo simulations as will be described in Sec. 2.2. The field autocorrelation function is calculated from g 1 ðτÞ ¼ hE Ã ðtÞEðt þ τÞi∕hjEðtÞj 2 i. For a semi-infinite medium, the analytical form of g 1 for CW-DCS system has been calculated as 7,9 Here, is the wavelength, k 0 is the central wave vector, C is the normalization constant, D B is the diffusion coefficient, and μ a and μ s 0 are the absorption and reduced scattering coefficients, respectively. In real time measurements, a simpler expression is preferred for fitting and Eq. (3) can be simplified to a single exponential decay at small τ limit 15,23 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 5 5 4 where τ c is the decay time constant that is utilized to quantify the tissue dynamics. Experimentally, instead of g 1 ðτÞ, the intensity correlation function g 2 ðτÞ ¼ hIðtÞIðt þ τÞi∕ hjIðtÞji 2 is measured, and g 2 ðτÞ is related to g 1 ðτÞ via the Siegert relation 24 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 4 8 7 with β ¼ 1 indicating complete coherence of the detected photons, and β < 1 accounting for loss of temporal coherence and/or detection of multiple modes such as the two polarization states of the electromagnetic field. The effect of the coherence length l c has been estimated as 25 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 4 0 7 Here, PðLÞ is the photon pathlength distribution. It is noteworthy that this relation is only valid when the incident light spectrum is a Gaussian function. For a semi-infinite highly scattering medium relevant for human brain measurements, PðLÞ can be calculated as 26,27 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 3 2 7 Here, v is the speed of light in the sample, τÞ is the pathlength-dependent normalized field temporal autocorrelation function E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 2 1 7 At τ ¼ 0, g 1 ðL; 0Þ ¼ 1, which provides the relation between β and l c in Eq. (6). The noise model for CW-DCS systems has been established analytically. 15 The standard deviation of ðg 2 ðτÞ − 1ÞÞ, σðτÞ, at each time delay τ is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 1 4 8 Here, T b is the bin time, m is the bin index, T is the measurement time window, hni is the average number of photons in a bin time For TD-DCS, a gating strategy is utilized to select a portion of the photons that arrive at the detector. For a gate delay time t s and gate width t w , only the photons that arrive within a measurement window are selected to calculate the electric field in Eq. (2), as shown in Fig. 1(b). Ideally, if both the widths of the IRF and the gate width are infinitely small, the gating strategy will only select photons with a single pathlength L ¼ t s Ã v, and g 1 ðτÞ is reduced to where v is the speed of the light in the tissue. However, both the IRF and the gate width will increase the width of the distribution of the photon pathlengths detected within a measurement window. The expression of g 2 ðτÞ for TD-DCS has been analytically obtained and the effect of the coherence length l c and IRF are taken into account. 18 At t s , the gate delay time dependent intensity autocorrelation function g 2 ðt s ; τÞ is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 6 1 6 Here the effective pathlength distribution Pðt s ; LÞ ¼ PðLÞI p ðt s − L∕vÞ is the pathlength distribution of photon trajectories measured at t s , where the normalized IRF profile is denoted as I p ðtÞ. The definition of coherence length is shown in Sec. 2.2, Step 4. It is noteworthy that the definition of the coherence length can be different in other studies. 18,25

Monte Carlo Wave Model
We have established a numerical recipe to simulate photon counts received at the detector as a function of time for both CW-DCS and TD-DCS measurements from first principles. The steps needed to construct the model are explained below.
Step 1: Simulate photon migration within a scattering medium using Monte-Carlo simulations.
The code for the Monte Carlo simulation is a derivation of that used in previous publications. 14,23,28 In the Monte Carlo simulation, we specify the source-detector separation and the optical properties of the sample including the scattering coefficient μ s and anisotropy factor g in different tissue types. We have used a two-layer sample geometry denoted as two tissue types in the Monte Carlo code as shown in Fig. 1(b). Here we have set g ¼ 0 and μ s ¼ μ s 0 ¼ 1 mm −1 for both tissue types in this study, but these can be easily configured to be different by researchers for future work. In biological tissue, the value of g obtained experimentally is typically above 0.9. 29 But in the diffusive regime as we are interested, the DCS results are only governed by the reduced scattering coefficient as indicated in Eq. (3). Thus, we can implement isotropic scattering with g ¼ 0 and use a smaller μ s by assigning μ s ¼ μ s 0 in the Monte Carlo simulation to improve computational efficiency, which is equivalent to g ¼ 0.9, μ s 0 ¼ μ s ð1 − gÞ. The effect of absorption will be taken into account at a later stage in Step 8. We will assign different blood flow dynamics to the two layers to demonstrate the DCS signal variation induced by blood flow change due to neural activation in the deeper layer. The number of scattering events N sni and pathlength L nti for the n'th photon in the i ti'th tissue type are recorded. The values of other parameters we have used in the model are specified in the caption of Fig. 1 Step 2: Calculate the arrival time t na of the n'th photon. We set the time at the peak of the IRF to be t ¼ 0. The transit time that a photon spent within the medium is t nL ¼ L n ∕v, where v ¼ c∕n is the speed of light within the medium, c ¼ 300 mm∕ns is the speed of the light in vacuum and the refractive index of the tissue is set to be n ¼ 1.33. The IRF time t nIRF is randomly drawn from the probability density distribution determined by the IRF profile, by obtaining the time at which a random number drawn within the interval [0 1] matches the cumulative probability. This takes into account the effect of the IRF on the distribution of the photon pathlengths detected within a specific measurement window. The arrival time of the photon is t na ¼ t nL þ t nIRF . In this paper, we have used an experimentally measured IRF as shown in Fig. 1(d). In the experiment, a PicoQuant VisIR laser at 765 nm was used as the source with a PPG512 to trim the pulse, and a PicoQuant PMA42 photomultiplier was used as the detector to reduce the long tail. Any IRF profile obtained experimentally can be utilized in this model to predict the performance for a particular TD-DCS system. Step 3: Select photons that fall within the specified measurement gate. For a gate delay time t s with a gate width t w , only the photons with t na that fall within the time span ½t s − t w ∕2; t s þ t w ∕2 are selected for the following calculations. It is noteworthy that if both the width of the IRF and the gate width t w are infinitely small, the photons with a single pathlength L n ¼ t s Ã v are selected, which is the ideal case described by Eq. (8), but is generally not achievable experimentally.
Step 4: For each photon, assign the spectrum of the input light that corresponds to a particular coherence length l c . Here we have utilized a Gaussian spectrum, and the relation between the Gaussian spectrum in real and k space is related via E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 6 0 6 Here, we have defined the coherence length l c as the width of the Gaussian. It is noteworthy that the definition of the coherence length can be different in other studies. We have combined the photon description of light with the wave description. In principle, we cannot assign a wavelength to a photon. However, each photon can be considered as a light trajectory for wave propagation, which allows us to calculate the accumulated phase along a light path for different wavelengths. In our model, any spectrum shape SðkÞ measured experimentally can be incorporated, which is one advantage of our numerical model over existing theoretical models, 18,25 where a Gaussian spectrum profile is always assumed.
Step 5: Assign an initial phase randomly drawn within ½0 2π for each photon trajectory ϕ n .
Here, Nð0;2D B ΔtÞ is a normal distribution with mean 0 and variance 2D B Δt. It is noteworthy that the motion of the scatterers or red blood cells is in general a combination of convective flow and shear induced diffusion. It has been experimentally observed and numerically demonstrated that the DCS signal is dominated by the diffusive behavior of RBCs red blood cells (RBCs). 14,30,31 Thus, we only model the diffusive motion of RBCs as in Eq. (12). The diffusion coefficient at the baseline, i.e. without brain activation, is set to be D B ¼ 1 Ã 10 −6 mm 2 ∕s in the model, which provides a decay time constant close to experimental measurements. 9,23 Step 7: Calculate the momentum transferq nis ¼ ðk out −k in Þ nis for each scattering event. Here i denotes the i'th wave vector k for the s'th scattering event of the n'th photon. The angle of a scattering event can be sampled from the phase function. However, since we are assuming g ¼ 0 where the scattering is isotropic as described in Step 1, we can simply generate unit vectorsñ out with a random direction and calculateñ out −ñ in using a coordinate space wherẽ n in ¼ ð0;0; 1Þ. The momentum transfer is thenq nis ¼ k i Ã ðñ out −ñ in Þ ns .
Step 8: At every time point t, obtain the contribution of the n'th photon and i'th k vector to the electric field as It is noteworthy that if μ a is different in different tissue types, μ a L n ¼ P ti μ ati L nti , where μ ati and L nti are the absorption coefficient and photon pathlength in the ti'th tissue type Step 9: Calculate the intensity fluctuations IðtÞ. This is calculated by summing over the contribution from all the photon trajectories and wavelengths E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 6 ; 5 9 1 IðtÞ ¼ Here we have ignored the term e iωt , where ω is the angular frequency of the light, and we have performed an incoherent sum over the k vectors since typically the measurement time is long enough such that light waves with different angular frequencies or wavelengths do not interfere. An example of the simulated intensity temporal fluctuations IðtÞ is shown in Fig. 2(a).
Step 10: Convert intensity IðtÞ to photon counts NðtÞ. If the average photon flux within a bin time Δt is N bin , we then normalize IðtÞ such that the average value of hIðt n Þi n ¼ N bin , where Iðt n Þ is the intensity within the n'th bin time t n . Then we perform a Poisson draw for each bin time as, Nðt n Þ ¼ PoisðIðt n ÞÞ. This process converts the continuous function IðtÞ to discrete values as well as incorporates shot noise into the numerically generated photon counts NðtÞ. Other sources of noise can also be added to NðtÞ, including dark noise, but this is out of the scope of this paper as dark noise is usually negligible in experimental systems. An example of the conversion from IðtÞ in Fig. 2(a) to NðtÞ is shown in Fig. 2(b). The Poisson draw here is not unique due to the nature of shot noise.
Step 11: Calculate the intensity autocorrelation function g 2 ðτÞ ¼ hNðtÞNðt þ τÞi∕hNðtÞi 2 . It is noteworthy that if we skip Step 10, IðtÞ can also be utilized to calculate the intrinsic noise in g 2 ðτÞ ¼ hIðtÞIðt þ τÞi∕hIðtÞi 2 that arises from speckle statistics due to a finite measurement time, 32 which does not include the effect of shot noise.
For CW-DCS simulations, we skip Steps 2 and 3, and use all the recorded photons to calculate the speckle fluctuations utilizing Eqs. (13) and (14). In Sec. 3, we compare our modeling results with existing theoretical predictions described in Sec. 2.1.

Results
We first simulate CW-DCS measurements with the measurement configuration described in Fig. 1(b). We compare our results with the theoretical predictions of both g 2 ðτÞ in Eq. (4) and the noise model in g 2 ðτÞ and σðτÞ in Eq. (9) as shown in Figs. 3(a) and 3(b). We see that our numerical results are in good agreements with the theoretical predictions. These results are obtained for a single wavelength at 800 nm under the assumption that the coherence length is infinitely long. Next, we simulate CW-DCS measurements using multiple wavelengths with central wavelength 800 nm to compare with the theoretical predictions of the reduction of β due to loss of coherence described by Eq. (6), by varying the width of the source spectrum (i.e., the coherence length). The pathlength distribution PðLÞ of all the photons detected at the detector agrees with the theoretical prediction in Eq. (7), as shown in Fig. 4(a). With the same PðLÞ, we have calculated β ¼ g 2 ð0Þ − 1 as a function of the coherence length l c using our numerical model and compared with the results calculated from the theoretical prediction in Eq. (6). In these calculations, we have used IðtÞ instead of NðtÞ to obtain g 2 ðτÞ in Step 11 of Sec. 2.2, since the effect of shot noise is not included in Eq. (6). 25 After validating our modeling results with existing theoretical predictions for CW-DCS systems, we now simulate TD-DCS measurements with and without the IRF profile in Fig. 1(d) taken into account, using a coherence length of l c ¼ 90 mm for the input light source. For the case without the effect of the IRF, the IRF profile is assumed to be infinitely narrow, i.e., t nIRF ¼ 0, t na ¼ t nL ¼ L∕v in Step 2 of Sec. 2.2. The gate delay time dependent pathlength distributions Pðt s ; LÞ for t s ¼ 0.5, 1.5, and 2.5 ns and g 2 ðτÞ with and without the IRF are shown in Figs. 5(a) and 5(b). We see that if the effect of the IRF is ignored, the pathlength distribution Pðt s ; LÞ is solely determined by the gates including gate delay t s and gate width t w . Compared to Fig. 5(a), the pathlength distributions in Fig. 5(b) are broadened and their shapes are determined by the IRF profile in Fig. 1(d). The corresponding results of g 2 ðτÞ are shown in Figs. 5(c) and 5(d). Compared to Fig. 5(c), the values of β ¼ g 2 ð0Þ − 1 in Fig. 5(d) are smaller due to a broader distribution of pathlengths, as we have expected. To validate the modeling results, we compare the g 2 ðτÞ curves obtained from the Monte Carlo wave model (solid lines) with those obtained from the analytical expression in Eq. (10) (dashed lines) as shown in Figs. 5(c) and 5(d). We see that the modeling and analytical results are in good agreement. It is noteworthy that for Figs. 5(a) and 5(c), we have assumed a finite coherence length with an infinitely small width of the IRF, which is generally not physically realizable. This numerical analysis serves as a sanity check against established concepts. Next, we perform a noise analysis when only the speckle noise is considered, i.e., in the high photon flux limit. We first obtain g 2 ðτÞ ¼ hIðtÞIðt þ τÞi∕hIðtÞi 2 and then obtain τ c from fitting using Eq. (4). The signal is obtained as the average of the τ c values, which is denoted as τ c for simplicity. Noise is obtained from the stdðτ c Þ by running 100 instances of the Monte Carlo wave model simulations using the baseline sample properties. The SNR is defined  Other parameters are the same as in Fig. 3.
as SNR ¼ τ c ∕stdðτ c Þ. To calculate the CNR for estimating changes in blood flow, we produce an activated state by changing the diffusion coefficient from the baseline value D B ¼ 1 × 10 −6 mm 2 ∕s to D B ¼ 2.25 × 10 −6 mm 2 ∕s in the second layer (Tissue 2) shown in Fig. 1(b). The value of τ c 0 , which is the average τ c at the activated states, is obtained from the average of 10 instances of the wave model. The contrast is Δτ c ¼ τ c − τ c 0 , and the CNR is CNR ¼ Δτ c ∕stdðτ c Þ. Here, we have used a large change of the diffusion coefficient to minimize errors in the contrast induced by noise to save computational time, since this paper is mainly focused on demonstrations of the model instead of calculating the realistic contrast due to brain activation. A smaller change of D B can be applied, while more configurations used for averaging are required to obtain both τ c and τ c 0 . The comparison of the performance of TD-DCS with and without the effect of IRF, and CW-DCS in the speckle noise limit are shown in Fig. 6. The value τ c decreases with increasing t s in Fig. 6(a) as expected, since photons with longer pathlengths that undergo more scattering events are selected at larger t s . The SNR increases with t s as shown in Fig. 6(b). This is because for a fixed measurement time (T ¼ 10 ms in this case), more speckle evolutions are averaged at longer t s because of the shorter τ c . To quantify the specificity to brain activation in the second layer in Fig. 1(b), we obtain the contrast with activation only in the second layer, i.e., the diffusion coefficient is increased from baseline value D B ¼ 1 × 10 −6 mm 2 ∕s to D B ¼ 2.25 × 10 −6 mm 2 ∕s in the second layer (Tissue 2), and compared with the case when brain activation is applied to the full medium, i.e., the diffusion coefficient is increased from baseline value D B ¼ 1 × 10 −6 mm 2 ∕s to D B ¼ 2.25 × 10 −6 mm 2 ∕s in the whole simulation region with its contrast denoted as Δτ c;h . We see in Fig. 6(c) that the specificity to brain activation in the second layer, i.e., Δτ c ∕Δτ c;h , increases with increasing t s , and that at large t s , Δτ c ∕Δτ c;h is larger for TD-DCS than for CW-DCS. This demonstrates the ability of TD-DCS measurements at larger t s to achieve better sensitivity to tissue dynamics at deeper layers. Thus, TD-DCS can tune the depth specificity by changing the measurement gate delays with a fixed source-detector geometry, which is convenient for technologies such as highdensity DCS measurements. Also, CNR is higher for TD-DCS compared to CW-DCS in this high-photon flux, speckle noise limit as shown in Fig. 6(d). Now, we consider the more realistic case for human brain measurements where the photon flux is very low and the DCS signal is shot noise limited instead of speckle noise limited. The photon flux difference between TD-DCS and CW-DCS measurements needs to be taken into account. We assume the incident photon flux to be 100; 000 photons∕s for CW-DCS at ρ ¼ 20 mm. We then scale the photon flux accordingly by the number of photons detected in a measurement gate in Step 3 for TD-DCS and the total number of the detected photons in Step 1 of the Monte Carlo simulation for CW-DCS, as shown in Table 1. For the low photon fluxes typical of human measurements, long measurement times T are needed to get estimates of g 2 ðτÞ with sufficient quality that can be fit to obtain an estimate of τ c . Running the wave model for long T is computationally expensive. Further, we want to keep T fixed to be 10 ms to facilitate comparison with prior examples. We thus utilize the following trick to improve computational efficiency that gives us an estimate of the shot noise contribution to stdðτ c Þ (i.e., stdðτ c Þ shot ) that is distinct from the finite speckle sampling noise contribution to stdðτ c Þ (i.e., stdðτ c Þ speckle ). We can then get the total noise from the sum of the variances. Our trick is to sample the wave model for T ¼ 10 ms to get IðtÞ. We then obtain N avg ¼ 1000 independent Poisson draws of this IðtÞ get 1000 instances of NðtÞ. We obtain g 2 ðτÞ for each instance of NðtÞ and average the 1000 correlation functions to get a more smooth g 2;avg ðτÞ. We then fit to get τ c;avg from g 2;avg ðτÞ. We repeat this process 100 times to obtain stdðτ c;avg Þ. Finally, we obtain stdðτ c Þ shot ¼ ffiffiffiffiffiffiffiffiffi N avg p Ã stdðτ c;avg Þ, since shot noise contribution is proportional to 1∕ ffiffiffiffiffiffiffiffiffi N avg p . With these photon flux, we then calculate the CNR for TD-DCS at ρ ¼ 10;20 mm with the effect of the IRF and coherence length taken into account, and CW-DCS at ρ ¼ 10;20; 30 mm as shown in Fig. 7(a). We see that for TD-DCS, the maximum CNR for the cases we have calculated occurs at ρ ¼ 10 mm, t s ¼ 1.5 ns, when t w is fixed to be 0.2 ns, but its CNR is still lower than CW-DCS at ρ ¼ 20 mm. The best performance for the cases we have simulated is CW-DCS at the largest source-detector separation of ρ ¼ 30 mm. Thus, when photon flux is taken into account, the relative performance of TD-DCS as compared to its CW-DCS counterpart has degraded. One way to increase the photon flux in TD-DCS is to increase the gate width t w . We see that by increasing t w , the CNR of TD-DCS at t w ¼ 0.5 ns has surpassed the performance of CW-DCS at ρ ¼ 20 mm. But at large t w , the specificity to a particular photon pathlength L and thus the ability to differentiate tissue depths decreases.

Discussion
We have developed a Monte Carlo wave model that is capable of simulating both TD-DCS and CW-DCS measurements from first principles. This model can be utilized to calculate the noise in TD-DCS and compare with CW-DCS measurements. We have validated our model with existing analytical expressions of g 2 ðτÞ for CW-DCS, TD-DCS, and the noise model of CW-DCS. Our numerical model is capable of simulating any spectral profile, while existing theories always assume a Gaussian spectrum. 18,25 Our system is capable of dealing with different brain activation geometries other than the two layer system we have specified in this manuscript. From the comparison of TD-DCS and CW-DCS, we have seen that the performance of TD-DCS is not necessarily better compared to CW-DCS, with the metric given by the CNR induced by tissue dynamics change in the deeper layer of a two layer sample, mainly due to the relatively low photon flux in TD-DCS since only a portion of the detected photons are selected. Methods that can increase the photon flux for a measurement gate, such as temporal focusing, could significantly improve the performance of TD-DCS, as we have seen that in the high photon flux limit, TD-DCS out-performs CW-DCS. The performance of TD-DCS is highly sensitive to the measurement parameters including gate delay time t s , gate width t w , and source-detector separation ρ. It also depends on the hardware properties including the wavelength and coherence length of the input laser light, the size of the detector, and the IRF profile. An analysis of the full parameter space is out of the scope of this manuscript. For a give TD-DCS measurement and brain activation geometry, the optimized parameters can be determined using the code of our Monte Carlo wave model. Our model will guide the future experimental design of novel TD-DCS systems.

Disclosures
Dr. Boas has consulted with Facebook Inc./Meta Platforms Inc. on topics related to DCS that ultimately led to the ideas presented in this paper. Dr. Boas's interests were reviewed and are managed by Boston University in accordance with their conflict of interest policies. Table 1 The photon flux used for all the simulation cases in Fig. 7.